module adv_mod

	use namelist_mod
	use mesh_mod
	use vert_coord_mod
	use time_mod
	! use process_mod
	use state_mod
	use tracer_mod
	use history_mod
	use time_mod, old=>old_time_idx, new=>new_time_idx
	use operators_mod
	use step_mod
	use operators_mod
	! 2d test
	use solid_rotation_test_mod     
	use deform_test_mod
	use stationary_vortex_test_mod
	use moving_vortex_test_mod
	! 3d test
	use deform_flow_test_mod
	use hadley_merid_test_mod
	use thin_cloud_test_mod

	implicit none

	interface 
		subroutine set_ic_interface(mesh, tracer)
			import mesh_type, tracer_type
			type(mesh_type), intent(in) :: mesh
			type(tracer_type), intent(inout) :: tracer
		end subroutine set_ic_interface

		subroutine get_exact_interface(mesh, tracer, time_in_seconds)
			import mesh_type, tracer_type
			type(mesh_type), intent(in) :: mesh
			type(tracer_type), intent(inout) :: tracer
		  real(8), intent(in) :: time_in_seconds
		end subroutine get_exact_interface

		subroutine set_wind_interface(mesh, state, time_in_seconds)
		  import mesh_type, state_type
		  type(mesh_type ), intent(in   ) :: mesh
		  type(state_type), intent(inout) :: state
		  real(8), intent(in) :: time_in_seconds
		end subroutine set_wind_interface

		subroutine step_adv_interface(old_state, new_state, old_tracer, new_tracer, time_in_seconds)
			import state_type, tracer_type
			type(state_type), intent(inout) :: old_state
			type(state_type), intent(inout) :: new_state
			type(tracer_type), intent(inout) :: old_tracer
			type(tracer_type), intent(inout) :: new_tracer
			real(8), intent(in) :: time_in_seconds
		end subroutine step_adv_interface
	end interface

	procedure(set_ic_interface   ), pointer :: set_ic    => null()
	procedure(get_exact_interface), pointer :: get_exact => null()
	procedure(set_wind_interface ), pointer :: set_wind  => null()
	procedure(step_adv_interface ), pointer :: step_adv  => null()

contains

	subroutine adv_init()

		integer i

		call mesh%init(lon_halo_width=40, lat_halo_width=2, lev_halo_width=2)
		call vert_coord_init(num_lev)
		call time_init()
		! call process_init()
		if (.not. allocated(state)) allocate(state(2))
		do i = 1, size(state)
			call state(i)%init(mesh)
		end do

		if (.not. allocated(tracer)) allocate(tracer(2))
		do i = 1, size(tracer)
		  call tracer(i)%init(mesh)
		end do

		call history_init()

		select case (test_case)
		! 2d test
		case ('solid_rotation')
			call solid_rotation_test_init()
			set_ic   => solid_rotation_test_set_ic
			set_wind => solid_rotation_test_set_uv
		case ('deform_case1')
			call deform_test_init(1)
			set_ic   => deform_test_set_ic
			set_wind => deform_case1_test_set_uv
		case ('deform_case2')
			call deform_test_init(2)
			set_ic   => deform_test_set_ic
			set_wind => deform_case2_test_set_uv
		case ('deform_case3')
			call deform_test_init(3)
			set_ic   => deform_test_set_ic
			set_wind => deform_case3_test_set_uv
		case ('deform_case4')
			call deform_test_init(4)
			set_ic   => deform_test_set_ic
			set_wind => deform_case4_test_set_uv
		case ('stationary_vortex')
			call stationary_vortex_test_init()
			set_ic    => stationary_vortex_test_set_ic
			set_wind  => stationary_vortex_test_set_uv
			get_exact => stationary_vortex_test_get_exact
		case ('moving_vortex')
			call moving_vortex_test_init()
			set_ic    => moving_vortex_test_set_ic
			set_wind  => moving_vortex_test_set_uv
			get_exact => moving_vortex_test_get_exact
		! 3d test 
		case ('dcmip11')
			call deform_flow_test_init()
			set_ic   => deform_flow_test_set_ic
			set_wind => deform_flow_test_set_wind
		case ('dcmip12')
			call hadley_merid_test_init()
			set_ic   => hadley_merid_test_set_ic
			set_wind => hadley_merid_test_set_wind
		case ('dcmip13')
			call thin_cloud_test_init()
			set_ic   => thin_cloud_test_set_ic
			set_wind => thin_cloud_test_set_wind
		case default
			print*, "Unknown test case " // trim(test_case) // '!'
			stop
		end select

		call step_init()

		select case (trim(advection_scheme))
		case ('split')
			step_adv => step_adv_split
		case ('nosplit')
			step_adv => step_adv_nosplit
		case ('ieva')
			step_adv => step_adv_ieva
		case default
			print*, "Unknown advection_scheme"
	  	stop
		end select

	end subroutine adv_init

	subroutine adv_run()

		integer nt, io_step

		io_step = 1

		call set_ic(mesh, tracer(old))
		if (test_case == 'stationary_vortex' .or. test_case == 'moving_vortex') then
		  call get_exact(mesh, tracer(old), real(0, 8))
	  end if
		call set_wind(mesh, state(old), elapsed_seconds)
		call diagnose(tracer(old), state(old), dt_adv)
	  call history_output(state(old), tracer(old), io_step)

	  nt = int(run_days * 86400._r8 / dt_adv)

	  do while (time_step <= nt)
	  	call set_wind(mesh, state(new), elapsed_seconds + dt_adv)
	  	call step_adv(state(old), state(new), tracer(old), tracer(new), dt_adv)
	  	if (test_case == 'stationary_vortex' .or. test_case == 'moving_vortex') then
	  	  call get_exact(mesh, tracer(old), elapsed_seconds)
	    end if
	  	call time_advance(dt_adv)
	  	call diagnose(tracer(old), state(old), dt_adv)
	  	if (mod(elapsed_seconds, history_interval) == 0) then
	  		print*, 'time_step=', time_step, 'of', nt
	  		io_step = io_step + 1
	  	  call history_output(state(old), tracer(old), io_step)
	    end if
	  end do

	end subroutine adv_run

	subroutine adv_final() 

		call mesh%clear()
		call vert_coord_final()
		if (allocated(state))  deallocate(state)
		if (allocated(tracer)) deallocate(tracer)
		call step_final()

	end subroutine adv_final

	subroutine diagnose(tracer, state, dt)

		type(tracer_type), intent(inout) :: tracer
		type(state_type) , intent(in) :: state
		integer i, j, k, itracer
		real(r8) tm
		real(8) dt
		logical exist

		do itracer = 1, size(tracer%q,4)
			tm = 0.0_r8
			do k = mesh%full_lev_ibeg, mesh%full_lev_iend
				do j = mesh%full_lat_ibeg, mesh%full_lat_iend
					do i = mesh%full_lon_ibeg, mesh%full_lon_iend
						tm = tm + tracer%q(i,j,k,itracer) * mesh%area_cell(j) * mesh%full_dlev(k)
					end do
				end do
			end do 
			tracer%tm(itracer) = tm
		end do

		call calc_hori_cfl(tracer, dt)
		call calc_vert_cfl(tracer%we, state%m_lev, tracer%cflz, dt)
		
		inquire(file="./max_cr_num.txt", exist=exist)
		if (exist) then
			open(11, file='./max_cr_num.txt', status="old", position="append", action="write")
		else
			open(11, file="./max_cr_num.txt", status="new", action="write")
			write(11, *) "cr_num_x    crnum_y    cr_num_z"
		end if
		write(11, *) maxval(tracer%cflx), maxval(tracer%cfly), maxval(tracer%cflz)
		close(11)

	end subroutine diagnose

end module adv_mod